
#include <libcmat.h>
#include <vector>
#include "gallop_common.h"
#include "2bwt_format.h"
#include "check_result.h"

#define CHECK_QUAL
#define CHECK_POS
#define CHECK_READ
#define CHECK_HIT
#define CHECK_STRAND


void soap2bin(AlignResult* alignResult, const unsigned int qid, 
		Soap_format* soap, const int readLength) {
	char* qscore = alignResult->get_qscore(qid);

	for(unsigned int i = 0; i < readLength; i++) {
		qscore[i] = soap->get_qual(i);
	}

	alignResult->_pos[qid] = soap->get_pos();
	alignResult->_hit[qid] = soap->get_hit();
	alignResult->_strand[qid] = soap->strand;
}

void readSoap(AlignResult* alignResult, ifstream& soapFile, const int readLength) {
	Soap_format soap;
	soapFile.clear();
	
	/*1. get the numAlign and readLength*/
	unsigned long numAlign = 0;
	for(std::string line; getline(soapFile, line);) {
		std::istringstream s(line);
		if((s>>soap) && (soap.get_read_len() == readLength)) {
			numAlign += 1;
		}
	}
	
	/*2. allocate the memory*/
	alignResult->clear();
	alignResult->alloc(numAlign, readLength);

	/*3. parse into the data structure*/
	soapFile.clear();
	soapFile.seekg(0);
	unsigned int count = 0;
	char* qscore = NULL;
	for(std::string line; getline(soapFile, line);) {
		std::istringstream s(line);
		if((s>>soap) && (soap.get_read_len() == readLength)) {
			soap2bin(alignResult, count, &soap, readLength);
			count++;
		}
	}
}

template<class T>
bool pairCompare(pair<unsigned int, T> a, pair<unsigned int, T> b) {
	if(a.first < b.first) {
		return true;
	} else if (a.first > b.first) {
		return false;
	} else {
		if(a.second < b.second) {
			return true;
		} else {
			return false;
		}
	}
}

template<class T>
void checkAttr(vector<pair<unsigned int, T> > &gold, vector<pair<unsigned int, T> > ref, string attr) {

        cout << "checking " << attr << " ..." << endl;
        stable_sort(gold.begin(), gold.end(), pairCompare<T>);
        stable_sort(ref.begin(), ref.end(), pairCompare<T>);
        bool pass = true;
        cout << "start ..." << endl;
	assert(gold.size() == ref.size());
        for(int i = 0; i < gold.size(); i++) {
                if(gold[i].second != ref[i].second) {
                        cout << "pos: " << gold[i].first << ", " << ref[i].first << endl;
                        cout << "gold: " << gold[i].second << endl;
                        cout << "ref:  " << ref[i].second << endl;
                        cout << endl;
                        pass = false;
			//break;
                }
        }
        cout << "***" << attr << " " << ((pass) ? "pass" : "failed") << "***" << endl;
}

bool compareSoapBin(const char* soapFileName, const char* galignerFileName, 
			const char* refFileName, const int readLength, 
			const char* dbsnpFileName = NULL) {
	bool isPass = true;

	/*read the binary alignment result*/
	AlignResultRead* read = new AlignResultRead(galignerFileName);
	cout << "open the GAligner binary output file done." << endl;

	/*open the reference and set it for AlignResultRead*/
	Sequence* ref = new Sequence;
	ref->read(refFileName);
	read->setRef(ref->getBase());
	cout << "chr_name: " << read->getChrName() << endl;	

	/*prepare for the sorted soap file*/
	ifstream soapFile(soapFileName, ifstream::in);
	Soap_format soap;

	/*prepare for the counting matrix generation*/
	//CMatrix* cmat = new CMatrix(refFileName, dbsnpFileName);

	/*compare the result*/
	cout << "start read and check the records of galigner ..." << endl;
	char* base = new char[readLength];
	char* qual = new char[readLength];
	char* readBase = new char[readLength];
	int hit, readLen;
	unsigned int pos;
	char strand;
	uchar2 offset; 
	char2 sub;
	unsigned int count = 1;

	//container for checking
	vector<pair<unsigned int, string> > qualGold, qualRef;
	vector<pair<unsigned int, string> > readGold, readRef;
	vector<pair<unsigned int, char> > strandGold, strandRef;
	vector<pair<unsigned int, int> > hitGold, hitRef;
	vector<pair<unsigned int, uchar2> > offsetRef;
	vector<pair<unsigned int, char2> > subRef;

	for(std::string line; getline(soapFile, line);) {
		//read soap
		std::istringstream s(line);
		if((s>>soap) && (soap.get_read_len() == readLength)) {
		} else {
			cout << "!!!ERROR: parsing soap file errors!" << endl;
			exit(-1);
		}


		//read binary file and check the result
		if(read->fetchOrdered(readBase, qual, &hit, &readLen, 
				&pos, &strand, &offset, &sub)) {
			
			if(readLen != soap.get_read_len()) {
				cout << "!!!! read length error" << endl;
				exit(-1);
			}

#ifdef CHECK_POS
			if(pos != soap.get_pos()) {
				cout << "!!! pos error: " << pos << ", " << soap.get_pos() << endl;
				isPass = false;
				exit(-1);
			}
#endif


#ifdef CHECK_READ
			/*for(int k = 0; k < readLength; k++) {
				if(readBase[k] != soap.get_base(k)) {
					cout << "!!! base error" << endl;
					cout << "line " << count << endl;
					cout << readBase << "\n" << soap.read << endl;
					if(soap.get_mismatch() == 1) {
						cout << soap.mis1 << endl;
						printf("%d, %c\n", offset.x, sub.x);
					} else if(soap.get_mismatch() == 2) {
						cout << soap.mis1 << ", " << soap.mis2 << endl;
						printf("%d, %c\n", offset.x, sub.x);
						printf("%d, %c\n", offset.y, sub.y);
					}
					cout << endl;
					isPass = false;
					break;
				}
			}*/
			readGold.push_back(pair<unsigned int, string>(pos, soap.read));
			readRef.push_back(pair<unsigned int, string>(pos, readBase));
#endif 

#ifdef CHECK_QUAL
			
			/*for(int k = 0; k < readLength; k++) {
				if(qual[k] != soap.get_qual(k)) {
					cout << "!!! qual error" << endl;
					isPass = false;
					cout << "position:  " << pos << endl;
					cout << "correct:   " << soap.qual << endl;
					cout << "incorrect: " << qual << endl;
			
					break;
				}
			}*/
			qualGold.push_back(pair<unsigned int, string>(pos, soap.qual));
			qualRef.push_back(pair<unsigned int, string>(pos, qual));

#endif


			if(read->getChrName() != soap.get_chr_name()) {
				cout << "!!!chr_name error" << endl;
				isPass = false;
				break;
			}

#ifdef CHECK_HIT
			/*if(hit != soap.get_hit()) {
				cout << "!!! hit error" << endl;
				isPass = false;
			}*/
			hitGold.push_back(pair<unsigned int, int>(pos, soap.get_hit()));
			hitRef.push_back(pair<unsigned int, int>(pos, hit));
#endif

#ifdef CHECK_STRAND
			/*if(strand != soap.strand) {
				cout << count << ": " << endl;
				cout << "!!! strand error: " << endl;
				printf("%c, %c\n", strand, soap.strand);
				isPass = false;
			}*/
			strandGold.push_back(pair<unsigned int, char>(pos, soap.strand));
			strandRef.push_back(pair<unsigned int, char>(pos, strand));
#endif
		} else {
			cout << "!!!!ERROR: parsing binary alignment file error" << endl;
			isPass = false;
		}

		//now we update the counting matrix
		//cmat->update(readLen, readBase, qual, pos, hit, strand, (read->getChrName()).c_str());

		count++;
	}


#ifdef CHECK_READ	
	/*for(int i = 0; i < readGold.size(); i++) {
		if(readGold[i].first == 137970) {
			cout << "gold: " << readGold[i].first << ": " << readGold[i].second << endl;
		}
		if(readRef[i].first == 137970) {
			cout << "ref:  " << readRef[i].first << ": " << readRef[i].second << endl;
		}
	}*/
	checkAttr<string>(readGold, readRef, "read");
#endif

#ifdef CHECK_QUAL
	checkAttr<string>(qualGold, qualRef, "qual");
#endif

#ifdef CHECK_STRAND
	checkAttr<char>(strandGold, strandRef, "strand");
#endif

#ifdef CHECK_HIT
	checkAttr<int>(hitGold, hitRef, "hit");
#endif

	//output the counting matrix
	//const char* cmatFileName = "/lhome/mian/data/chr21.test.cmat";
	//cmat->write(cmatFileName);

	//delete gold;
	delete read;
	//delete cmat;
	soapFile.close();
	ref->clear();
	
	return isPass;
}

void compareSoapBin_api(int argc, char** argv) {
	if(argc != 5 && argc != 6) {
		cout << "1. the name of **SORTED** soap output" << endl;
		cout << "2. the name of galigner output" << endl;
		cout << "3. the name of the reference" << endl;
		cout << "4. the read length" << endl;
		cout << "5. the dbsnp file name (optional)" << endl;
		exit(-1);
	}

	const char* soapFileName = argv[1];
	const char* galignerFileName = argv[2];
	const char* refFileName = argv[3];
	const int readLength = atoi(argv[4]);
	char* dbsnpFileName = NULL;
	if(argc == 6) {
		dbsnpFileName = argv[5];
	}

	if(compareSoapBin(soapFileName, galignerFileName, 
			refFileName, readLength, dbsnpFileName)) {
		cout << "*** Test passed ***" << endl;
	} else {
		cout << "*** Test FAILED !!! ***" << endl;
	}
}

